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Abstract 

The Leonids show meteor storms in a period of 33 years, and known as one of the most active meteor 
showers. It has recently shown a meteor stream consisting of several narrow dust trails made by meteoroids 
ejected from a parent comet. Hence, an analysis of the temporal behavior of the meteor flux is important 
to study the structure of the trails. However, statistical inference for the count data is not an easy task, 
because of its Poisson characteristics. We carried out a wide-field video observation of the Leonid meteor 
storm in 2001. We formulated a state-of-the-art statistical analysis, which is called a self-organizing state 
space model, to infer the true behavior of the dust density of the trails properly from the meteor count 
data. From this analysis, we found that the trails have a fairly smooth spatial structure, with small 
and dense clumps that cause a temporal burst of meteor flux. We also proved that the time behavior 
(trend) of the fluxes of bright meteors and that of faint meteors are significantly different. In addition we 
comment on some other application of the self-organizing state-space model in fields related to astronomy 
and astrophysics. 

Key words: comets: individual (comet 55P/Tempel-Tuttle) — interplanetary medium — meteors, 
meteoroids — methods: statistical 



1. Introduction 

1.1. Leonid Meteor Storm as a Stochastic Counting 
Process 

The Leonids show meteor storms in a period of 33 years, 
and known as one of the most active meteor showers. It 
is believed that the Leonid meteor storm occurs when the 
Earth approaches the orbit of the dust trails of the comet 
55P/Tempel-Tuttle. Many authors have predicted that 
the Leonid meteor storm would occur in 2001 (see e.g., 
Watanabe et al. 2002 and references therein) . It has been 
recognized that a meteor stream consists of several narrow 
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dust trails, each of which is made by meteoroids ejected at 
a particular return of the parent comet (e.g., Watanabe et 
al. 1997; Brown et al. 2002). Hence, the temporal behavior 
of the meteor flux provides important information of the 
spatial and density structure of the dust trails. 

We carried out a wide-field video observation of the 
Leonid meteor storm in 2001 (Shiki et al. 2003), and found 
the following results through a statistical analysis: 

1. The time variation of the hourly rate (HR) shows 
many small peaks. 

2. The HR profile of bright meteors and that of faint 
meteors are significantly different. 

3. The small peaks arc associated with a burst of faint 
meteor flux. 
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4. Some signatures of Poisson-like process are found 
in the interval distribution of meteors and in their 
autocorrelation function. 

Many previous studies tried to analyze the time variation 
of the meteor flux (e.g., Watanabe et al. 1997 and refer- 
ences therein). However, it is not a trivial task to infer 
the true behavior of the trail density, since it suffers from 
a strong statistical fluctuation. The difficulty is due to 
the fact that measuring the meteor flux is a kind of typ- 
ical counting process, whose fluctuation is basically not 
Gaussian, but is characterized by a Poisson distribution. 
The Poisson counting process has a variance equal to its 
mean value, i.e., when the mean is large, the associating 
variance is also large. This character prevents us from ap- 
plying the popular classical time series analysis methods 
to the count data, because some basic assumptions are 
utterly violated (e.g., Brockwcll, Davis 1996; Cameron, 
Trivedi 1998). 

1.2. Development of Time Series Analysis 

Classical time series analyses have often appeared in 
an astronomical context. Data from astronomy as well 
as from other physical, biological, or econometric studies 
consist of a sequence of numbers, {xi,X2,X3, • • • ,xt}, ob- 
tained by measuring the quantity xt during a sequence of 
times. The subscripts represent discrete time that runs 
from t = 1 to T. This discrete treatment of time is suf- 
ficient for a variety of practical applications. 1 A compre- 
hensive summary of the classical time series analysis and 
its applications to astronomical datasets can be found in 
Scargle (1981). 

Today, in order to handle a wider range of time series, 
including some latent variables, the state space model has 
been proposed from the field of engineering and system 
optimal control. Especially, the merit of using the state 
space form is that it can properly treat time- varying pa- 
rameters in the system. We note that classical time se- 
ries models can be defined as special cases of the gen- 
eral state space model. The well-known Kalman filter 
(Kalman 1960) gives an algorithm to estimate the system 
parameters recursively, i.e., it gives one-step-ahead esti- 
mations (called 'filtering') everytime we have a new data 
point in the series. 

A Kalman filter assumes the Gaussianity for the system 
noise terms, which is not suitable for the data we often 
encounter in astronomical applications. Further, as seen 
in section 3.1, the Kalman filter is designed for a linear 
system. A generalized state space model with nonlinear- 
ity and non-Gaussianity has been proposed to overcome 
these shortcomings (Kitagawa, Gersch 1996, and refer- 
ences therein). An excellent guide for their applications 
can be found in Brockwell, Davis (1996). 

However, there still remains an additional difficulty 
concerning astronomical data analysis. Astronomers fre- 
quently meet a situation in which they must handle a 

1 Here, we concentrate on a set of equally sampled time series 
data, but, as mentioned later, the method presented in this ar- 
ticle can also be applied to irregularly sampled data. 



dataset with small counts, such as a very low-level signal 
of the faintest sources. Count data introduce complica- 
tions of discreteness and heteroskedasticity. 2 The inclu- 
sion of zero counts appears to be a pitfall to apply, e.g., 
usual regression methods. Despite its frequency that we 
should tackle such datasets, time series models for count- 
ing data are in their infancy, yet remarkably many models 
have been developed. Cameron, Trivedi (1998) provides a 
through discussion on general counting problems, includ- 
ing time series counting data. 

Now we close the chronicle of time series analysis. 
The development is concisely summarized in Kitagawa, 
Sato (2001). In this paper, we propose a suitable method 
of analyzing time series count data, which sometimes have 
zeros in the sequence. This approach, the self-organizing 
state space model, has been developed only very recently 
(Kitagawa 1998; Higuchi 1999, and references therein). It 
makes extensive use of a large computing power of modern 
computers. 

In this paper we statistically formulate the count data of 
the Leonid meteor storm. The count rate obeys a Poisson 
process with a latent, time- varying Poisson intensity, At. 
We made an attempt to estimate the temporal behav- 
ior of the hidden parameter A* from the observed count 
data. We should note that our data have some gaps in ob- 
servations caused by the length of video tapes and other 
reasons. Our method can easily overcome such gaps, and 
properly infer the value in the observational time gaps. 

The rest of the present paper is as follows. In section 3, 
we formulate a self-organizing state space model of the 
Leonid count data. We start with a linear state space 
model, and then develop toward more general methods. 
We present the results and discussions in section 4. sec- 
tion 5 is devoted to a summary. 

2. Data 

The original data consist of the magnitudes and ob- 
served time of the meteors, which were recorded on a video 
tape. Although they include some sporadic meteors and 
meteors belonging to other meteor showers, here we con- 
centrate only on the Leonids in the analysis. Detailed 
data descriptions are found in Shiki et al. (2003). 

In the usual manner of the meteor-shower analysis, an 
hourly rate (HR) is used to describe the time variation of 
the meteor flux. However, in this study, it is more ap- 
propriate to present the data based on the count rate per 
minute; hence, we use the count rate [min -1 ] throughout 
this paper. 

3. Method 

3.1. State Space Model 

A classical state space model consists of the following 
equations: 

x t = F t x t -i + Gtvt , v t ~N(Q,Q t ), (1) 
y t = H t x t + w t , wt~N(p,Rt). (2) 

2 This technical term means that the variance is not constant. 
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Here, Xt is called a state vector, which represents the (un- 
observed) state of the system, and y t is the observed se- 
quence of data. The idea underlying the model is that the 
development of the system over time is determined by Xt 
according to equation (1). However, because x t cannot 
be observed directly, we must base the analysis on obser- 
vations, y t . We call equation (1) the state equation, and 
equation (2) the observation equation. The error terms, Vt 
and Wt, are distributed according to Gaussian probability 
distribution functions, N(0,Q t ) and iV(0,Rt), respectively, 
where Qt and R t are covariance matrices. 

Matrices F t , Gt, H t , Q t , and Rt are initially assumed to 
be known, and the error terms, Vt and Wt, are assumed to 
be serially independent and independent of each other at 
all time. In practice, some or all of the matrices depend 
on the elements of an unknown parameter vector, 9. 

A considerable advantage of the state space approach 
is the ease with which missing observations can be dealt 
with. The estimation problems in time series analysis can 
be classified into the following three categories with re- 
spect to the dependence on the observed data {yi, ■■■ ,yt} 
to estimate the state vector xt'. 

1- Vi:t-i = {yi,---,yt-i} ■ prediction, 

2. yi-.t = {yi,---,yt} ■ filtering, and 

3- yi:u = {Vi,---,Vv}{u>t) : smoothing. 

Kalman recursion equations give a one-step-ahead pre- 
diction of the state vector, xt, and its error covariance 
matrix, by using {yi, ■ ■ ■ ,y t _i}, and filter the series when 
we have new data, yt- We do not go any further into the 
details of the Kalman filter here. For implementation, see, 
e.g., Harvey (1981), Brockwell, Davis (1996) and Durbin, 
Koopman (2001). If we have a missing in data sequence, 
we simply perform prediction without filtering, and go 
to the next step. This point is extensively discussed by 
Akaike, Kitagawa (1999). 

3.2. Generalized State Space Model 

We can consider a nonlinear non-Gaussian state space 
model as being an extension of the linear case (Kitagawa 
1987): 



x t = F t (x t -i,vt) 
y t = H t (x t ,w t ) . 



(3) 
(4) 



Again, the first is the state equation and the second is 
the observation equation. These times, v t and Wt, are 
the system and observation noise with non-Gaussian den- 
sities, qt(vt) and r t (w t ), respectively. The initial state 
Xo is assumed to be distributed with the probability den- 
sity Po(xo). Functions F t {x,v) and H t (x,w) are nonlinear 
ones of the state vector and noise (Brockwell, Davis 1996; 
Durbin, Koopman 2001; Kitagawa, Sato 2001). 

It is convenient to express the model in a general form 
based on the conditional distributions: 



Xi 

Vi 



Qt{-\xt-i, 
Rt(-\x t ) . 



(5) 
(6) 



For general state space models, the conditional distribu- 
tions become non-Gaussian and their distributions cannot 
be completely specified by the mean vectors and the co- 
variance matrices, which is different from the case of a 
Gaussian linear state space model and a Kalman filter. A 
non-Gaussian filter is expressed as follows: 



p(xt\yi-.t-i) = J p(x t \xt-i)p(xt-i\yi:t-i)dx t -i , (7) 

, I \ p{yt\xt)p{xt\yi:t-\) , a , 

p{xt\yi:t) = 1—\ ^ ' ( 8 

p{y t \yi:t-i) 

where p(yt\yi-.t-i) is the predictive distribution of yt, 

p(yt\yi-.t-i) = / p(y t \xt)p(xt\yi:t-i)dx t . (9) 



However, a direct implementation of the formula requires 
computationally intense numerical integration that can 
only be feasible for some limited case. 

3.3. Monte Carlo Filter 

Instead of approximating the distribution and per- 
forming heavy numerical integration, we can use Monte 
Carlo filtering, i.e., producing a large number of real- 
izations which can be extracted from the distribution 
(Kitagawa 1996). This method needs much less computa- 
tional power. 

The procedure is as follows: 

1. Generate a random number, Xq ^pq(xq), for j = 
1,---,N. 

2. Repeat the following steps for t = 1, • • • , T: 

(a) Generate a random number v\ ~ q(v) for j = 
l,-,N. 

(b) Compute p[ 3) = F[x { ^\ , u t 0) ] for j = 1, ■ ■ ■ ,N. 

(c) Compute w^' = p[y t \p[^] for j = 1, • • • ,N. 

(d) Generate x[^ (j = 1, • • • , N) by resampling of 
pP (j = l,---,N). 

In step 2(b), = F[x[_i,Vt ] can be considered to be 
independent realizations from the predictive distribution, 
p{xt\y\:t-i) ■ Given the observation y t and the realization 
p*f\ we obtain the importance weight, u>j [step 2(c)], 
i.e., the likelihood with respect to the observation, y t . The 
posterior probability of the realization is 

Prob[x t =p[ j) |y 1:t ] 

= Prob[a: t =pr\yi:t-i,Vt] 

p[y t \p[ j) ]P r ob[xt=pi j) \yi:t-i} 



'Zk=iP[yt\P ( t K '] p r° } °[xt =Pt K> \yi-.t-i] 



(fe)i 



(j) 



EN ( 
fc=l W t 



(10) 



This means that the cumulative distribution function of 
Prob[:ct =p[^\yi-t] can be expressed by a step function, 



With this general state space model, we can handle 
discrete- valued time series as well as discrete-state models. 
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Prob[x t =p < t ) \yi:t] = 



N 



E 



N 

k=l 



W 



Y w t 'I[x,pf'] , (11) 3.5. Model for the Leonid Meteor Storm 



iPt 



X 1 ) 



XN) 



with step sizes w\ u>i 

which has a unit 



which has jumps at pi , • 
Here, I(x, z) is the Heaviside function, 
jump at x = z. 

For the next step of the prediction, we must repre- 
sent this distribution function by an empirical distribution 
with the form 



Prob[x t 



-Yi 



N 



N 



x,Pi 



This can be done by resampling of {p 
probabilities, 



Prob[xi 



■■P { t j) \yi:t] 



E 



N 
k=l 



W 



(/■■) 



(12) 



Pt } with 



(13) 



We calculate this step using the von Neumann's 
acception-rejection method (see e.g., Knuth 1998 for de- 
tails). 

One of the main purposes of the time series analysis is 
to estimate some hidden parameters, 9, from the observed 
data. The likelihood of the time series model specified by 
the parameter 9 is obtained by 



L (8) = p(yi,---,yr) = IJp(j/t|2/i:i-i; 



(14) 



Here, we formulate our problem to estimate the tempo- 
ral trend of the Leonid meteor shower hidden by a statisti- 
cal fluctuation, which may be Poissonian. For the Poisson 
count process, the observation equation is expressed as 



y t ~ Poisson(A t ) 



-At 



t 



T 



(Higuchi 1999; Higuchi 2001; Kitagawa, Sato 2001). 
system model becomes quite simple, as 



.2 



logCT 2 



(19) 



The 



(20) 



where \Xt '■ 
such that 



log At. We adopt a first-order smooth trend 



+ v t , v t ~N(0,o-l). 



This treatment 
enables us to handle an arbitrary trend, because it only 
assumes that the first-order difference of the trend is 
small. The parameter a 2 (logcr 2 ) is simultaneously es- 
timated by a recursive procedure. 3 As recommended by 
Higuchi (2001), we set log a\ ~ {/([-6.0,-2.0]) as the ini- 
tial distribution of a 2 ^ , where U([a,b]) denotes the uniform 
distribution between a and b. 

As already mentioned above, Leonid meteors seem to 
behave like a Poisson process locally, and hence this model 
is appropriate to describe this data (see Shiki et al. 2003). 

4. Results and Discussions 



Along with the above discussion, we can approximate the 
conditional density by 



p{yt\y\-.t- 



N 

N ^ 



,(./) 



(15) 



3.4- Self- Organizing State Space Model 

In principle, the maximum likelihood estimate of 9 
is obtained by maximizing the log-likelihood, log L(6). 
However, in practice the sampling error and long com- 
putational time often renders the direct maximum like- 
lihood method impractical. To remedy this prob- 
lem, Kitagawa (1998) proposed a sophisticated method. 
Instead of estimating the parameter 9 by the maximum 
likelihood, we consider a Bayesian estimation by augment- 
ing the state vector, x t , with an unknown parameter, 9, 
as 



Xi 



(16) 



and construct the state space model for z t as 

z t = F*(z t . 1 ,v t ), (17) 

Vt~R{-\zt). (18) 

This is called 'the self-organizing state space model'. In 
equation (16), the parameter is unknown, but constant 
such that 9 t — 9 t -i = ■ ■ ■ = 9. This can immediately ex- 
tended to the time- varying (hyper)parameter case, but we 
do not go any further here. 



The observations were made from 14 h 41 m to 20 h 03 ra UT 
(322 minutes), hence T = 322 in equation (19). 

4--1- Global Behavior 

We present the count rate data of the Leonid meteor 
storm and its estimate for the true density distribution in 
Figure 1. At a glance we can see that the density esti- 
mate has a fairly smooth spatial structure, and that the 
violent statistical fluctuation is significantly suppressed by 
the self-organizing state space method. The estimation 
uncertainty, which is simultaneously estimated by the re- 
cursive estimation process, is ~0.1. The ±l-er uncertainty 
envelopes are also plotted in Figure 1, but it is hard to 
resolve on the figure. This result shows that most of the 
'spikes' in the meteor count data are merely a consequence 
of the Poisson fluctuation, and have no physical substance. 
Therefore, even if a spike appears to be strong, it may be 
explained by a large variance (standard deviation) of the 
Poisson process. 

By our method, only those spikes which cannot be re- 
garded as a mere fluctuation are detected in the trend 
estimate. The most prominent feature is the burst of me- 
teor flux just before the first global peak of the storm. We 
can observe some other small peaks at 100, 130, 160, 210, 
250, 260, and 270 min in Figure 1. They clearly corre- 
spond to the peaks suggested by an analysis of Shiki et 
al. (2003). Interestingly, the first spike of the count at 



The logarithm is defined as logx = log 10 x. It is set to conserve 
the positive definiteness of the trend. 




Fig. 1. Count rate data of the Leonid meteor storm and its 
estimate for the true density distribution. The vertical dotted 
lines indicate the boundaries of the data. Data arc missing 
in small intervals between the dotted lines. The observed 
counts are shown by the thin histogram, and the thick zigged 
curve presents an estimate of the true count behavior, At. 
On the right-hand side we show the HR (hourly rate) for the 
convenience of readers. 



Fig. 2. Count rate data of the Leonid meteor storm for me- 
teors brighter than 3 mag, and its estimate for the true den- 
sity distribution. Again, the vertical dotted lines indicate the 
boundaries of the data. The observed counts are shown by 
the thin histogram, and the thick zigged curve presents an 
estimate of the true count. 



45 min (~ 15 25 m ) appeared to be a sudden increase of 
the underlying meteor flux, and not a spiky burst of count 
in the estimate. 

Thus, we conclude that the true temporal trend of the 
dust trail which caused the Leonid meteor storm is glob- 
ally smooth, with small and dense clumps associated with 
bursts of meteor counts. This confirms the suggestion 
from a classical analysis by Shiki et al. (2003), and pro- 
vides a statistically rigorous basis on it. 

4-2. Magnitude Dependence 

Next, we divided the data into two classes, bright and 
faint samples, and applied the self-organizing state space 
method to both of them, just as we did for the whole sam- 
ple. We set the boundary of the bright and faint samples 
at 3 mag. The results of the bright and faint samples are 
shown in Figures 2 and 3, respectively 

A drastic difference is found between the counts of the 
bright and faint samples. The most striking feature is 
that there is no trend of bursts in bright sample counts at 
200 min (~ 6 h 00 m UT), whereas a prominent burst exists 
in the faint counts. Other weaker bursts also stem from 
the faint count behavior. In other words, the bright sam- 
ple count is relatively smooth and its variation is small, 
while the faint count temporarily varies with a very sim- 
ilar trend of the total count profile. This clearly shows 
that most of the bursts have been dominated by meteors 
fainter than 3 mag. 

We also find a clear excess of bright meteors to faint 
ones after 270 min (19 h 20 m UT). It is an unexpected 
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Fig. 3. Same as Figure 2, except that it is for meteors fainter 
than 3 mag. 
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Fig. 4. Zenith-corrected meteor flux estimates. 

trend, because it is widely accepted that fainter meteors 
are more numerous than brighter ones. One may suspect 
the effect of the elevation of the radiation point, but this 
result is unchanged by the correction because the correc- 
tion has no magnitude dependence. This suggests a bias 
in the distributions of larger and smaller meteorites, which 
may be the origins of brighter and fainter meteors, respec- 
tively. 

Hughes (1973) reported a clear difference in the cumu- 
lative influxes of meteors brighter and fainter than 3 mag. 
Our result may be closely related to his conclusion, but 
we do not go further here. 

So far, we have not corrected the zenith effect on the 
meteor flux so as to avoid any unnecessary intricacy in 
statistical modeling. If one may wish to have a zenith- 
corrected count rate (or cquivalcntly, ZHR), wc should 
merely make a correction of the obtained estimate here. 
Figure 4 shows the corrected meteor flux and ZHR, to- 
gether with those of bright and faint subsamplcs. For the 
correction, wc simply multiplied 1 / sin h (h : elevation 
of the radiation point) to the flux estimates. If we use 
an empirical formula of the form 1/ (sin ft,) 7 (7 ~ 1.4: see 
e.g., Jenniskens 1994), the trend around 15 h -17 h would be 
more emphasized. We note that the sum of the estimates 
for these subsampcls perfectly agrees with the estimates 
for a whole sample. This also shows that our approach 
is very powerful, robust, and consistent for this type of 
analysis. In figure 4, the coincidence of the burst spikes 
in the whole sample and faint subsample is impressive. 

4.-3. Future Prospects of the Self- Organizing State Space 
Model for Astrophysical Applications 

Before closing this article wc would like to devote a 
subsection to some future prospects for applying the self- 
organizing state space model approach. This approach 



has a very wide range of its applicability: for example, it 
can be used to estimate an extremely faint optical source 
variability, and to analyze count rates of low-level pho- 
tons, X-ray or cosmic ray detectors, which are regarded 
as representatives of typical count processes. 

Another important aspect is its robustness against ir- 
regular sampling. Hence, we can easily apply this method 
to the photometric sequence data of gravitationally lensed 
objects to measure their time delay in variability. 

Higuchi (2001) illustrated an interesting application to 
an estimation of spiral density wave in Saturn's ring ob- 
served by Voyager, which is known to have a varying frcq- 
nency along with a radial position (Horn et al. 1996). He 
applied the self-organizing state space model approach to 
the data and showed a beautiful result. Sunspot num- 
ber data are also very popular count process in statistical 
science (Higuchi 1999). 

Thus, wc expect a variety of applications of this anal- 
ysis. Now, with this approach, we do not have to 
worry about unrealistic assumptions of stationarity nor 
Gaussianity that are hardly expected for real datasets, in 
spite of the fact that they are often required for the pop- 
ular classical time series analysis. Also, we will never be 
annoyed by irregularly appearing observational gaps, sam- 
pling inhomogeneity, or hetcroskedasticity that arc often 
inherent in various astronomical datasets. 

5. Summary and Conclusions 

The Leonids show meteor storms in a period of 33 years, 
and known as one of the most active meteor showers. It 
has recently shown a meteor stream consisting of several 
narrow dust trails made by meteoroids ejected from a par- 
ent comet; hence, an analysis of the temporal behavior of 
the meteor flux is important for studying the structure 
of the trails. However, statistical inference for the count 
data is not an easy task, because of its Poisson charac- 
teristics. We carried out a wide-field video observation of 
the Leonid meteor storm in 2001 (Shiki et al. 2003). 

In this study, we formulated a state-of-the-art statisti- 
cal analysis, which is called the self-organizing state space 
model, to infer the true behavior of the dust density of 
the trails properly from the meteor count data. From this 
analysis we found that the trails have a fairly smooth spa- 
tial structure, with small and dense clumps, which cause 
a temporal burst of meteor flux. We also confirmed that 
the time behavior (trend) of the fluxes of bright and faint 
meteors are significantly different. 

In summary, for the first time we obtained a reliable 
estimate of the true dust trail density profile by the self- 
organizing state space approach. 

First of all, we are greatly indebted to Prof. Tomoyuki 
Higuchi, the referee, whose careful reading and comments 
imporved the quality and rigor of this paper much. T.T.T. 
also thanks Prof. Peter Brockwcll for developing and pro- 
viding their analysis software ITSM, which has brought a 
number of insights into our pre-analysis, and Dr. Takako 
T. Ishii for detailed instruction of the coding for the devel- 
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opment of the self-organizing state space model. T.T.T. 
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